#ifndef __NURBS_PATCH_NURBS__
#define __NURBS_PATCH_NURBS__

#include "element.h"
#include "primitive.h"
#include "math/bsplines.h"

class SNurbs : public ElementImplementation<SNurbs> {
public:
	SNurbs(const boost::multi_array<Node*, 2>& c, const Knotvector& U, const Knotvector& V, size_t p, size_t q, Attribute& attribute, size_t id);
	
	virtual const Node* GetNode(size_t i) const;
	virtual Node* GetNode(size_t i);
	const Node* GetNode(size_t i, size_t j) const;
	Node* GetNode(size_t i, size_t j);

	template<typename T> void EvaluateImplementation(T& x, const Vector& P, size_t attributeHandle) const;
	template<typename T> void IntegrateImplementation(T& x, const Node* node, size_t attributeHandle) const;
	
	std::pair<Vec2n, Vec2n> Domain() const {return domain;}
	bool InDomain(const Vec2n& x) const;
	Vec3n Normal(const Vec2n& x, const std::string& attribute) const;

	virtual Matrix Stiffness(const std::string& attribute) const;
	virtual Matrix Mass(const std::string& attribute, const bool lumped) const;
	virtual Vector Force(const std::string& attribute) const;
	
	virtual size_t SizeNodes() const {return c.shape()[0] * c.shape()[1];}

	size_t GetDegreeU() const {return p;}
	size_t GetDegreeV() const {return p;}
	const Knotvector& GetKnotvectorU() const {return U;}
	const Knotvector& GetKnotvectorV() const {return V;}
	size_t GetSamplesU() const {return B[0][0][0].size();}
	size_t GetSamplesV() const {return B[1][0][0].size();}
	
	nfloat BasisFun(const math::Vec2ui& index, const math::Vec2ui& D, const math::Vec2ui& k) const;

	nfloat BasisFun(const size_t index, const size_t dim, const size_t degree, const size_t i) const {
		return B[dim][degree][i][index];
	}

	boost::array<nfloat, 2> Parameter(const math::Vec2ui& index) const {
		boost::array<nfloat, 2> param = {parameter[0][index[0]], parameter[1][index[1]]};
		return param;
	}

	nfloat Weight(const math::Vec2ui& index) const {
		return qweight[0][index[0]] * qweight[1][index[1]]; 
	}
	
	Matrix Jacobian(const Vec2n& param, const std::string& attribute, bool useWeights) const;
	Matrix Jacobian(const math::Vec2ui& index, const std::string& attribute, bool useWeights) const;
	
	bool Init();

	std::pair< std::vector<Vec3n>, std::vector<math::Vec3ui> > Tesselate(const size_t eppU, const size_t eppV, const std::string& attribute) const;
	std::pair< std::vector<Vec3n>, std::vector<math::Vec3ui> > Tesselate(const Array2V2n& param, const std::string& attribute) const;

	virtual void Dispatch(Dispatcher* d) {
		d->Visit(this);
	}

	bool IsExtraordinary(const Node* node) const {
		return node->Get<size_t>("extraordinary") != 0;
	}
	boost::multi_array<nfloat, 2> QuadraturePoints() const;
private:
	boost::multi_array<Node*, 2> c;
	Knotvector U, V;
	size_t p, q;
	std::pair<Vec2n, Vec2n> domain;
	mutable std::map<size_t, size_t> mapGlobalToLocal;
	boost::array< boost::tuples::tuple<bool, nfloat, nfloat>, 4> bc;

	nfloat Stiffness(const math::Vec2ui& k0, const math::Vec2ui& k1, const std::string& attribute, bool useWeights) const;
	nfloat Mass(const math::Vec2ui& k0, const math::Vec2ui& k1, const std::string& attribute, bool useWeights) const;
		
	nfloat BasisFun(const size_t i, const size_t j, const std::vector<nfloat>& bu, const std::vector<nfloat>& bv) const;
	nfloat BasisFun(const Vec2n& param, const math::Vec2ui& k) const;
	Vec2n Grad(const Vec2n& param, const math::Vec2ui& k, bool useWeights) const;
	Vec2n Grad(const math::Vec2ui& index, const math::Vec2ui& k, bool useWeights) const;
	
	Matrix Stiffness(const math::Vec2ui& index, const nfloat lambda, const nfloat mu, const bool nonlinear, const std::string& materialAttribute, const std::string& positionAttribute) const;
	Matrix Beta(const math::Vec2ui& index, const std::string& materialAttribute) const;
	Matrix Deformation(const math::Vec2ui& index, const std::string& materialAttribute, const std::string& positionAttribute) const;
	std::pair<Matrix, nfloat> Force(const math::Vec2ui& index, const nfloat lambda, const nfloat mu, const std::string& materialAttribute, const std::string& positionAttribute) const;
	nfloat Mass(const size_t i, const size_t j, const nfloat density, const std::string& materialAttribute) const;

	math::Vec2ui LocalIndex(const Node* v) const;
	
	std::vector< std::vector< std::vector< std::vector<nfloat> > > > B;
	boost::array< std::vector<nfloat>, 2> qweight;
	std::vector< std::vector<nfloat> > parameter;
	void EvaluateBasis();
};

template<typename T>
void SNurbs::EvaluateImplementation(T& x, const Vector& P, size_t attributeHandle) const {
	assert(P.size() == 2);
	
	std::vector<nfloat> bu(p + 1), bv(q + 1);
	for (size_t i = 0; i <= p; i++) bu[i] = math::spline::BasisFun<nfloat>(P(0), i, p, U);
	for (size_t i = 0; i <= q; i++) bv[i] = math::spline::BasisFun<nfloat>(P(1), i, q, V);

	x = T(0.0);
	if(Exist("weight")) {
		size_t weightHandle = AttributeHandle("weight");
		nfloat weight = 0.0;

		for (size_t i = 0; i <= p; i++) {
			for (size_t j = 0; j <= q; j++) {
				nfloat w = GetNode(i, j)->Get<nfloat>(weightHandle);

				x += w * GetNode(i, j)->Get<T>(attributeHandle) * BasisFun(i, j, bu, bv);
				weight += w * BasisFun(i, j, bu, bv);
			}
		}

		x = (1.0 / weight) * x;
	} else {
		for (size_t i = 0; i <= p; i++) {
			for (size_t j = 0; j <= q; j++) {
				x += GetNode(i, j)->Get<T>(attributeHandle) * BasisFun(i, j, bu, bv);
			}
		}
	}
}

template<typename T>
void SNurbs::IntegrateImplementation(T& x, const Node* node, size_t attributeHandle) const {
	math::Vec2ui index = LocalIndex(static_cast<const Node*>(node));

	boost::multi_array<nfloat, 2> gqp = QuadraturePoints();
	const size_t n = gqp.shape()[0];

	bool useWeights = Exist("weight");

	x = T(0.0);
	for(size_t i = 0; i < n; i++) {
		for(size_t j = 0; j < n; j++) {
			const Vec2n param = math::V2n(gqp[i][0], gqp[j][1]);
			const nfloat& wi = gqp[i][2];
			const nfloat& wj = gqp[j][2];
		
			const Matrix J = Jacobian(param, "position", useWeights);
			const Vec3n u = math::V3n(J(0, 0), J(1, 0), J(2, 0));
			const Vec3n v = math::V3n(J(0, 1), J(1, 1), J(2, 1));
			const nfloat Jdet = length(cross(u, v));

			Vec3n P = Evaluate<Vec3n>(param, "position");
	
			std::vector<nfloat> bu(p + 1), bv(q + 1);
			for (size_t k = 0; k <= p; k++) bu[k] = math::spline::BasisFun<nfloat>(param[0], k, p, U);
			for (size_t k = 0; k <= q; k++) bv[k] = math::spline::BasisFun<nfloat>(param[1], k, q, V);

			x += (wi * wj) * Evaluate<T>(param, attributeHandle) * BasisFun(param, index) * Jdet;
		}
	}
}

#endif //__NURBS_PATCH_NURBS__
